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Preface 

I A  f  -s  '/A*'5;;  is 

'V 

This  study  investigated  the  performance  of  an  on-board  filter  algorithm  in 
estimating  the  relative  positions  of  each  element  in  a  k*n  satellite  ^cluster.  The 
estimator  chosen  was  the  U-D  covariance  factorization  version  of  the  Kalman  filter. 

The  performance  was  evaluated  by  comparing  the  square-root  of  the  position 
covariance  eigenvalues  to  the  magnitude  of  the  true  position  errors.  Even  though 
filter  tuning  was  limited  to  the  diagonals  of  the  dynamics  noise  matrix,  the  position 
errors  were  well  within  the  accuracy  requirements  for  all  test  cases. 

In  programing  and  writing  my  thesis,  I  received  invaluable  help  from  my 
faculty  advisor,  Dr.  William  Wiesel.  His  patience  and  guidance  are  greatly  appre¬ 
ciated.  I  also  thank  Dr.  Peter  Majbeck  for  his  counsel.  Special  thanks  to  my  wife 
Susan  and  son  David  for  their  understanding  and  support  over  the  past  months. 
Finally,  I  wish  to  thank  the  Lord  who  gave  me  the  strength  to  keep  going  during 
the  long  hours. 

Michael  L.  P.  Ward 
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Abstract 


AFIT/GA/AA/88D-13 

The  relativex  position  determination  of  the  satellites  in  an  orbital  cluster  is 
investigated.  Thetm- board  estimator  is  the  U-D  covariance  factorisation  ^filter 
with  dynamics  based  on  the  Clohessy- Wiltshire  equations.  Performance  is  mea¬ 
sured  by  comparing  the  square  roc .  of  the  position  covariance  eigenvalues  to  the 
magnitude  of  true  position  errors.  True  errors  are  also  compared  to  minimum  accu¬ 
racy  requirements.  Filter  tuning  is  limited  to  adjusting  the  diagonal  entries  of  the 
dynamics  noise  matrix.  Test  cases  include:  perfect  initial  conditions  with  perfect 
range  data;  and  perfect  initial  conditions  with  noise  corrupted  measurements. 
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ESTIMATED  SATELLITE  CLUSTER  ELEMENTS 
IN  NEAR  CIRCULAR  ORBIT 


\b-  I.  Introduction 

A 

The  arraying  of^  satellites  into  clusters  has  been  under  consideration  for  a' 
number  of  years  for  commercial  communications  applications.  Stationing  these 
satellites  in  a  restricted  section  of  an  orbit  allows  for  in-orbit  spares  to  replace 
satellites  as  they  fail,  thereby  dramatically  increasing  overall  system  reliability. 
Reliability  is,  of  course,  of  major  concern  for  military  applications  as  well.  An 
additional  benefit  to  ib^military  planners  increased  wartime  system  survivability.  4 
A  cluster  basing  strategy  allows  for  the  use  of  less  expensive,  less  sophisticated 
satellites.  Survivability  is  enhanced  because  the  loss  of  a  satellite  does  not  destroy 
the  system. 

One  ‘proposed  use  of  such  a  satellite  cluster  would  be  a  space  based  radar. 

In  order  to  form  a  cohesive  image  of  the  radar  return,  the  relative  positions  of  the 
satellites  need  to  be  known  to  within  one  quarter  of  the  radar  wavelength(  [3:page 
443].  For  the  purposes  of  this  thesis,  25  meters, will  be'  considered  sufficiently 
accurate, (see  appendix  A). 

A  number  of  authors  have  published  studies  in  the  area  of  satellite  clusters. 
The  idea  was  first  introduced  by  P.  S.  Visher,  who  was  granted  a  patent  in  1983 
[2:page  1713].  Many  authors  have  done  follow-up  studies  on  cluster  geometries 
and  station  keeping  [6,2].  Few  have  addressed  the  stochastic  nature  of  the  position 
determination  problem. 

One  writer  who  has  addressed  this  problem  is  John  Murdoch  [8].  Specifi¬ 
cally,  he  examined  the  possibility  of  on-board  relative  position  determination  of 


the  satellites.  However,  he  performed  a  least-squares  analysts  which  is  less  appro¬ 
priate  than  a  recursive  filter  for  on-board  applications.  Also,  his  cluster  consisted 
of  only  two  satellites  in  near-geostationary  orbits.  An  actual  cluster  would  very 
likely  contain  many  more  satellites. 

This  thesis  is  not  concerned  with  a  specific  cluster  geometry  nor  with  station 
keeping  but  with  the  study  of  a  recursive  filter  for  on-board  relative  position  de¬ 
termination.  The  cluster  consists  of  up  to  ten  satellites  in  near-circular,  low-earth 
orbit.  To  complete  the  study,  a  truth  model  based  on  the  two- body  equations 
of  motion  was  develo,  <*d  to  provide  measurement  data  to  a  simulated  on-board 
Kalman  Filter.  A  single  sample  Monte  Carlo  analysis  wu  performed  to  rate  filter 
performance  by  comparing  the  covariance  standard  deviation  to  the  true  error. 


II.  Problem  Description 


t.l  Introduction 

In  the  satellite  cluster  concept,  each  satellite  would  estimate  the  positions  of 
every  element  of  the  array  to  get  its  own  position.  A  pulse  used  to  synchronise 
the  on-board  clocks  of  the  satellites  could  be  used  to  provide  range  data  x  for  the 
estimator.  For  the  remainder  of  this  study,  the  filter  algorithm  will  be  designed  for 
satellite  one  but  csm  be  generalised  for  the  other  satellites.  Figure  2.1  illustrates 
the  basic  concept  and  will  be  referred  to  in  later  analysis. 

In  order  to  test  the  concept  of  an  on-board  filter,  it  was  necessary  to  build 
a  real  world  or  truth  model  agaiust  which  the  filter  performance  could  be  judged. 
Since  the  purpose  of  this  thesis  was  to  evaluate  the  concept,  the  truth  model  was 
simplified  to  the  two-body  equations  of  motion.  The  truth  model  outputs  the 
true  state  x,(t)  and  the  measurements  *t(f).  These  measurements  are  used  by  the 
satellite  Kalman  filter  (figure  2.2)  to  output  an  estimated  state  x(t).  The  critical 
elements  of  the  filter  states  y(t)  are  then  compared  to  the  corresponding  elements 
of  the  true  state  y(t)  after  measurement  update  to  arrive  at  the  true  error  e,(t). 

e<  =  y  -  yt  (2.1) 

This  true  error  can  then  be  compared  to  the  square  root  of  the  appropriate  eigen¬ 
values  of  the  covariance  matrix  P  to  see  if  the  filter  performs  as  well  as  it  believes 
it  is  performing  [4:page  339]. 

t.l  Truth  Model 

The  truth  model  is  based  on  a  group  of  up  to  ten  satellite?  orbiting  the  earth 
as  a  (ingle  unit.  In  order  to  maintain  cluster  integrity,  they  have  to  remain  within 
a  prescribed  spatial  volume.  This  means  all  satellite  orbital  periods  must  be  equal 
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Figure  2.2.  Performance  Evaluation  of  a  Kalman  Filter 


or  the  formation  will  drift  apart.  Therefore,  the  model  first  creates  the  cluster  and 
then  propagates  each  satellite  through  its  orbit. 

As  stated  earlier,  the  output  of  the  truth  model  is  relative  range  data,  for 
input  to  the  Kalman  filter,  and  true  position,  for  comparison  to  the  estimated 
position.  Since  the  measurement  data  is  relative  information,  the  absolute  position 
cannot  be  estimated.  Instead,  position  is  described  with  respect  to  a  rotating 
reference  point  as  shown  in  Figure  2.1 

This  point  is  the  origin  of  a  reference  frame  (x,y,z)  rotating  with  respect  to 
an  earth  centered  inertial  frame  (t,j,fc)  at  a  constant  angular  rate  u>  k.  Note  (i,j) 
and  (x,y)  are  always  coplanar  and  that  k  is  always  parallel  to  z.  The  reference 
point  is  in  a  circular  orbit  at  radius  R  and  at  circular  orbit  velocity  of 


Vc,  = 


where  y  is  the  gravitational  parameter  of  the  earth.  The  initial  position  of  the 


reference  point  in  the  fix- A  frame  is 


rr«/ 


R  0  0 


(2.2) 


with  an  initial  velocity  of 

Vr  */ 


(2.3) 


Since  this  thesis  is  not  concerned  with  a  specific  geometry,  the  satellites  of  the 
cluster  are  positioned  at  random  about  the  reference  point  but  within  a  specified 
radius  rad : 


r,  = 


R 

0 

0 


+  n<  rad 


(2.4) 


where  i  =  1, 2, . . . ,  a  (the  total  number  of  satellites  in  the  cluster)  and  where  n<  is 
a  particular  vector  realization  of  a  uniform  random  number  generator  with  outputs 
between  -0.5  and  0.5.  Appropriate  values  for  the  velocity  vector  v,  can  be  found 
from  the  Clohessy- Wiltshire  equations  [9]  (these  equations  will  be  explained  in 
detail  in  the  next  section): 


Vi  i  =  rj(n  i  -  rrtf  i) 
Vi  k  =  rjri  k 


(2.5) 


v’here  q  is  the  mean  motion.  As  previously  stated,  the  velocity  must  be  specified 
so  that  the  period  of  all  satellites  are  equal.  Since  the  period  depends  only  on  the 
size  of  the  semi-major  axis  a  [l:page  33],  this  criterion  can  be  used: 

,2 
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E 

Ti 


ZE 

2a 


(2.6) 


Solving  for  Vi  j  yields 


j  =  (2^(~  ~  ^-)  ~  vi  *  -  v]  fc)> 


(2.7) 
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with  the  positive  square  root  chosen  so  that  each  satellite  will  maintain  formation 
integrity. 

By  using  this  position  and  velocity  information,  the  input  for  the  Kalman 
filter  can  be  formed.  Let  xt(<  =  0)  be  the  initial  true  state  defined  as  the  position 
and  velocity  of  each  satellite  in  the  rotating  frame  (ROT).  Initially,  the  axes  of 
both  coordinate  frames  are  aligned  therefore  the  position  of  each  satellite  in  the 
rotating  frame  is 

r<[ROT]  =  (r«  ~  r*-e/)[F/X]  (2-8) 

and  the  velocity  is 

vqnor]  =  vi[fia]  -  “>  x  r*[FJX)  (2.9) 

where  oj  =  ve,/R.  The  velocity  can  also  be  expressed  in  terms  of  relative  position 
as 

v«[rot]  =  (v*  -  vr</)jFJjrj  -  u>  x  r^oT]  (2.10) 

therefore  the  state  at  t  =  0  is 


xe(0)  = 


ri[Ror] 

vt[Ror] 


r«[flo:r] 

v,[Ror] 


(2.11) 


The  other  desired  output  is  a  j  —  1  dimensional  measurement  vector  of  the 

form: 


it 


ri  -  r2| 


|ri  -  r. 


•+■  ut 


(2.12) 


where  ut  is  a  zero-mean  white  Gaussian  noise  with  covariance  Rt  [4:page  330].  This 
noise  is  the  best  guess  of  the  real  world  errors  in  the  measurement  process.  The 
range  measurements  are  computed  from  impulse  signals  sent  to  each  satellite  by 
every  other  cluster  element  at  specified  times.  The  major  errors  result  from  clock 
synuironization  differences  between  satellites  as  well  as  the  resolution  limitations 
of  this  on-board  clock  [8:page  1016]. 

For  this  study,  noise  sources  axe  not  specified  individual1  y  but  are  lumped 
together  into  a  total  generic  measurement  error  err.  This  err  is  multiplied  by  the 
output  of  a  Gaussian,  random  number  generator  with  statistics  of  mean  of  0  and 
standard  deviation  of  ±1.0.  Thus,  Rt  is  an  (s  —  l)-by-(s  —  1)  matrix  with  diagonal 
entries  of  err2  and  off-diagonal  entries  of  zero.  Note  the  off-diagonal  zeros  imply 
that  the  measurements  are  independent;  this  assumption  is  unrealistic  given  this 
problem  statement  but  adequate  for  a  proof-of-concept  study.  Some  other  possible 
error  sources  are  biases  caused  by  misalignments  in  tracking  sensors.  These  biases 
are  not  addressed  in  this  thesis. 

The  next  task  of  the  truth  model  is  determining  the  future  positions  and 
velocities  of  the  satellites.  This  requires  solving  the  Kepler  problem  using  the  / 
and  g  equations  in  terms  of  eccentric  anomaly  E  [l:page  219].  This  solution  begins 
with  determining  additional  initial  conditions. 

Assume  the  initial  position  r0  and  velocity  v0  of  each  satellite  is  known  in 
the  inertial  frame.  Then  the  initial  eccentric  anomaly  Eo  and  mean  anomaly  Mo 
for  each  satellite  can  be  found  by  the  following  procedure: 


=  u  [(W®  ~~  TM>)  r#  (r°  *  Vo)  V® 


(2.13) 


where  e  is  the  eccentricity  vector  (constant  for  each  satellite) 


r-% 
\  ■• 


where  v  i  s  the  true  anomaly 


En  =  cos 


_!  f  (  +  COS  Up  \ 
\  1  +  «  COS  Vp  / 


and  finally' 


Mo  —  En  —  €  sin  En 


(2.15) 


(2.16) 


With  these  initial  conditions  in  hand,  the  Kepler  problem  can  now  be  solved. 
The  Kepler  problem  is  simply  the  solution  to  the  equation  [l:page  220] 

M  =  E  —  <  sin  E 

at  some  future  time  t.  M  at  this  time  is  easily  found  by 

M  —  \!  — ~  t  —  2  kir  +  Mo 


(2.17) 

r 

where  A;  is  an  integer  number  equal  to  the  number  of  times  the  satellite  has  passed 
through  M0.  Since  E  cannot  be  solved  for  directly,  it  must  be  determined  through 
an  iteration  scheme.  The  technique  chosen  for  this  thesis  was  a  Newton  iteration 
method.  Starting  with  a  first,  guess  of  E„  =  M,  the  corresponding  value  of  M„ 
using  Eq  (2.16)  is 

Mn  =  En  —  t  sin  En 


resulting  in  a  new  guess  (En+x)  of 


En+l  —  En  + 


M  - 


(2.18) 


1  —  €  COS  En 

This  iteration  continues  until  the  difference  M  —  Mn  becomes  sufficiently  small  and 
thus  determines  the  final  value  of  the  eccentric  anomaly  Ef. 

Now,  the  r  and  v  vectors  of  each  satellite  can  be  found  [l:page  198]: 


r  =  /  r0  +  g  v0 
v  =  /  r0  4-  g  v0 


(2.19) 

(2.20) 
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where  /,5,/,5  are 


/  =  1 - (1  —  cos  A E) 

r0 


(A E  -  sin  A E) 


.  y/fTa  Bin  AE 
r  r0 

g  =  1  —  —  (1  —  cos  A  E) 


(2.21) 


(2.22) 


(2.23) 

(2.24) 


The  two  unknown  values  in  the  above  equations  (AE  and  r)  are  determined  by 


A  E  =  Ef  —  Eq 
r  =  o(l  —  ecos  Ef) 


(2.25) 

(2.26) 


Since  r  and  v  of  each  satellite  are  now  known,  the  relative  position  to  the 
reference  point  in  the  rotating  frame  can  be  found.  First,  the  current  location  of 
the  reference  point  in  the  inertial  frame  can  readily  be  determined  by: 


Jicos0 

rre/(f7*]  =  iZsintf 


(2.27) 


— vet  sin  8 


Vrc f{FIX)  =  Vc,COS  9 


(2.28) 


where  9  =  u>t.  The  relative  position  and  velocity  is  then 


rrei [F/JT]  =  (r(t)  -  I*re/(03  =  rl*  +  r2j  +  ^ 


(2.29) 


Vrel[FIX]  =  [v(<)  -  Vre/(t)}  =  Vii  +  v2j  +  v3k 


(2.30) 
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Figure  2.3.  Cluster  State 

Unlike  the  situation  in  Eqs  2.8  and  2.10,  the  inertial  and  rotating  coordinate  frames 
are  not  aligned.  Therefore,  a  coordinate  transformation  must  be  performed  result¬ 
ing  in  the  following  equation: 


rq  cos  0  +  r j  sin  0 
—rx  sin  0  4-  r2  cos  0 
?z 

Vi  cos  0  4-  fj  sin  0  4-  u>  r2 
— v\  sin  0  +  Vj  cos  0  —  uri 


r[ROT] 

V[nor] 


(2.31) 


The  true  state  vector  is  formed  according  to  Eq  2.11  and  illustrated  in  Figure  2.3. 
The  out  of  plane  position  z  is  not  shown  but  is  measured  in  the  same  way  as 
the  illustrated  x  and  y  components.  Finally,  the  measurement  vector  is  formed 
according  to  Eq  2.12. 


2-9 


t.S  Kalman  Filter 

Now  that  the  ‘real  world’  model  has  been  completed  to  provide  the  nec¬ 
essary  data  for  the  filter,  the  filter  itself  can  be  developed.  Assume  the  system 
description  is  in  the  form  of  a  linear  stochastic  differential  equation  (to  be  devel¬ 
oped  later)  which  describes  state  propagation;  with  discrete-time,  noise  corrupted, 
nonlinear  measurements  as  available  inputs.  The  appropriate  estimator  is  the 
extended  Kalman  Filter.  Extended,  because  the  nonlinear  measurements  are  lin¬ 
earised  about  the  new  state  estimate  at  each  sample  time  as  opposed  to  being 
linearized  about  a  fixed  nominal  state  for  all  time  [5:page  42]. 

Since  the  dynamics  are  linear,  the  standard  Kalman  filter  propagation  equa¬ 
tions  can  be  used  [4:page  220]. 

*(<r)  =  *((.'.  (2-32) 

P(tr)  =  +GJ((i_,)qJ(i,.1)Gj(li.1)  (2.33) 

where 

Qd  The  covariance  of  the  dynamics  driving  noise. 

x(t; )  The  estimated  state  after  time  propagation. 

x(<+)  The  estimated  state  aft:-  the  measurement  update. 

Gd  Equals  the  identity  matrix  (because  the  model  is  an  equivalent  discrete¬ 
time  representation  of  a  continuous-time  system  [4:page  377]). 

$  The  state  transition  matrix. 

For  convenience,  the  time  argument  will  be  dropped  for  the  remainder  of 
this  discussion  mid  (tf)  will  be  written  as  (±).  The  first  time  propagation  oc¬ 
curs  prior  to  any  measurement  updates,  therefore  the  initial  conditions  x0  and  P0 
must  be  determined.  As  a  first  cut,  let  Xo  equal  xt(0)  from  the  truth  model  with 
the  position  measured  in  kilometers  and  velocity  in  kilometers  per  second.  Also, 
let  Po  be  diagonal  with  position  and  velocity  covariances  of  order  10'JJbmJ  and 
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10 ~,(km/sec)i  respectively.  These  values  result  from  an  assumed  initial  position 
accuracy  of  10  m  and  a  corresponding  velocity  error  of  2  m/sec.  This  Po  can  be 
replaced  with  the  steady-state  value  once  the  filter  has  been  tested.  Let  equal 
the  squared  difference  between  Xt  and  x  after  the  first  time  propagation.  A  sample 
time  of  350  seconds  was  chosen  to  ensure  the  relative  satellite  motion  could  be 
observed. 

Because  of  the  nonlinearity  of  the  measurements,  the  extended  Kalman  filter 
update  equations  must  be  used  [5:page  44]. 


K  =  P(-)Hr{  HP(-)Hr  +  R /  }“* 

(2.34) 

x(+)  =  x(-)  +  K{z-h[*(-)]} 

(2.35) 

P(+)  =  P(-)-KHP(-) 

(2.36) 

K  is  the  KaT  a  filter  gain  and  R/  is  the  filter  measurement  noise  covariance 
equaling  R(  from  the  truth  model. 

For  an  understanding  of  the  h  vector,  refer  back  to  Eq  2.12.  The  expectation 
of  this  equation  can  be  written  as 

*  =  h(x(<)) 

Therefor*,  h(x,f)  is  the  filter’s  estimate  of  range  and 

z  -  h[x(  — )] 

is  the  residual  or  true  error  in  range  used  to  update  the  state.  Thus,  h  has  the 
form 

v/(*  1  -  *a)J  +  (l/i  -  l/j)J  +  (*i  -  2j)s  hi 

h  =  ;  =  ;  (2.37) 

\/{ri  -  ®.)J  +  (yi  -  y.)3  +  (zi  -  zt)7  ht. i 
The  matrix  K  comes  from  linearizing  the  li  vector  and  evaluating  it  at  x(  — ): 
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Which  has  the  form: 


“ss 


Xaft(-) 


(2.38) 


H  = 


Hx 

Ht 

% 

ff.-i 


-Hx 

0 


0 


0 

0 


(2.39) 


where 


Hx  = 


vizjh  ss.  0  0  0 


(2.40) 


#,_1  =  >£=*1.  0  0  0 

1  fc*-i  *a-t 


(2.41) 


2.3.1  First  Filter.  The  first  filter  designed  Bind  tested  was  the  extended 
Kalman  filter  described  above.  The  only  element  in  the  algorithm  that  has  not 
been  fully  explained  in  this  text  is  the  ♦  matrix.  The  state  transition  matrix  was 
developed  from  the  Cloliessy- Wiltshire  equations  of  motion  [9:page  3]  as 


*  —  2rjy  —  3rjix  =  0 

ji  +  2i7*  =  0 
Z  +  TJ*Z  =  0 


(2.42) 

(2.43) 

(2.44) 


These  equations  can  be  integrated  about  the  initial  conditions  xo,xo>J/o»  etc. 
at  t  =  0  to  obtain  the  position  and  velocity  solutions  [9:pages  3-5] 


x(f) 


=  —  ^~yo  + 


*o  2 

yo  +  3xo  )  cos  rjt  4-  —  sin  rjt  +  4x0  — j/o 


(2.45) 
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y(t)  =  yo  -  (3ys  +  617*0)  t  +  I  —  +  6*0  I  sin  q <4-  —  cos  q< - (2.46) 

V  »?  /  n  n 


x(<)  —  to  count  +  ~  sin  f># 
tj 

*(t )  =  (2yo  +  317*0)  sin  qt  +  x0  cos  rjt 
y(t)  =  -3jfo  -  6q*o  +  (6q*o  -i-  4yo)cosqt  -  2*0  sin  qt 
i(<)  =  -x0qainq<  +  i0  cos  qi 


(2.47) 

(2.48) 

(2.49) 

(2.50) 


Notice  the  initial  downrange  position  yo  only  appears  as  an  additive  constant  in 
Eq  2.46  and  no  where  else.  This  will  be  important  in  later  discussion.  Finally,  the 
state  transition  matrix  can  be  formed  [9:page  6]. 


4  —  3  cos  0 

6(sin  V>  —  V1)  1 

0  0 

3»7  sin  V1  0 

6t7(cos  V’  —  1)  0 


J(1  -co«V>) 


COS  V1 


— r/sin 


l(coui>-l)  *sinifr-^  0 

0  0 

n 

cos  2  sin  V1  0 

— 2sinV’  —  3 +  4  cos  0 

0  0  cos  V’ 


(2.51) 


where  t/’  is  q6<  and  is  the  sample  time.  Using  this  $  matrix,  the  state  propagation 
equation  for  satellite  one  can  be  written  as 


Xi(-)  =  *iX!(  +  ) 


(2.52) 


Generalising  this  equation  for  3  satellites  gives 


Not?  tlte  ♦  matrix  ia  block  diagonal. 

One  of  the  possible  problems  associated  with  the  standard  Kalman  filter 
that  has  been  illustrated  thus  far  is  numerical  instability.  The  truth  model  was 
implemented  on  the  computer  in  double  precision  resulting  in  numerical  accuracy 
to  sixteen  significant  figures.  Since  the  filter  was  intended  to  operate  on  board  a 
simulated  satellite,  it  was  constrained  to  single  precision  or  eight  significant  figures. 
A  method  for  analysing  the  stability  of  a  given  filter  is  through  use  of  a  condition 
number  k( A)  [4:399].  The  condition  number  of  a  matrix  A  is 

*(A)  =  —  (2.54) 

in 

where  and  are  the  maximum  and  minimum  eigenvalues  of  Ar  A.  One 
of  the  indicators  of  instability  is  a  covariance  matrix  with  negative  eigenvalues 
because,  by  definition,  the  covariance  is  positive  semidefinite  [4:page  90]. 

When  this  filter  was  actually  implemented,  the  covariance  matrix  did  in  fact 
have  negative  eigenvalues.  Upon  investigation  of  k(  A),  it  was  discovered  that  the 
minimum  and  maximum  eigenvalues  of  PrP  were  on  the  order  of  10~lT  and  10* 
respectively.  This  resulted  in  a  condition  number  on  the  order  of  10*.  Studies  show 
numerical  difficulties  will  be  experienced  as  the  condition  number  approaches  10*, 
where  N  is  the  number  of  significant  digits.  It  is  evident  that,  as  described,  this 
filter  cannot  be  implemented  in  single  precision. 

t.S.t  Second  Filter.  One  method  of  solving  the  numerical  problems  is  using 
the  ‘U-D  covariance  factorisation'  version  of  the  Kalman  filter  [4:pages  392-399]. 
This  form  of  filter  achieves  twice  the  precision  with  the  same  wordlength  [4:400]. 
This  performance  results  from  factoring  the  covariance  into 


P(-)  =  U(-)D(-)Ur(-) 

(2.55) 

P(+)  =  U(+)D(+)Ur(  +  ) 

(2.56) 

where  the  U  matrix  i»  unitary  upper  triangular  anti  the  D  matrix  is  diagonal. 
This  method  guarantees  the  positive  semidefiniteness  of  the  covariance  and  the 
numerical  stability  and  accuracy  of  the  filter. 

The  algorithm  begins  with  the  same  initial  conditions  that  were  used  in  the 
Kalman  filter.  The  first  step  is  to  reduce  the  n-by-n  P0  (where  n  is  the  number 
of  states,  in  this  case  n  =  6s)  into  the  factored  form  UDUr  through  the  following 
steps:  First,  for  the  nth  column 


Dnn  —  Inn 
1 


uin  = 


t  =  n 

P,n/Dnn  i  =  n-  l,n-2,...,l 


(2.57) 


Then  for  the  remaining  columns,  j  =  n  —  l,n-2,  ...,1 


Djj 


Ua  = 


p»  -  z:.^Dkkujk 
0  i>  j 

1  *  —  J 

[■Pij  —  Sfcsj+i  DhhUikUjii^  / Djj  i  =  j  —  1,  j  —  2, . . . ,  1 


(2.58) 


This  procedure  can  still  be  used  if  the  Qj  and  R /  matrices  are  non-diagonal  (to 
be  discussed  later). 

Since  the  initial  values  for  U  and  D  have  been  determined,  the  state  can  be 
propagated  to  the  first  update  time.  Begin  with  setting  U(-f )  =  Uo  and  D(  +  )  = 
Do.  Then  form  the  n-by-2n  matrix  Y(  — )  as 

Y(-)  =  (  $U(  +  )  |  Gj  ]  (2.59) 

where  $  and  Gj  are  as  previously  defined.  Finally,  form  the  2n-by-2n  matrix 
D(-)by 


f  D(  +  )  0  1 


/  n  \ 


The  propagation  begins  by  forming 


Yr(-) 


ai  a  2 


(2.61) 


and  then  calculating  the  following  relationships  for  k  =  n,n  — 


C* 

=  D(  - )afc  ( ckj  =  Djj(  - )akJ ,  j  =  1 , 2, . . . ,  2n) 

Aek(-) 

T 

=  c* 

—  c  k/Dkk(-) 

(2.62) 

ujk(~) 

=  ajdfc  j  =  1,2,...,*:  -  1 

<  )®fc  J  —  1)  2, . . . ,  k  1 

The  final  portion  of  the  U  —  D  covariance  factorization  filter  algorithm  is  the 
scaler  measurement  update.  Let  the  first  row  of  Eq  2.39  be  the  1-by-n  matrix  Hi 
and  let  zx  be  the  corresponding  scaler  measurement  of  Eq  2.12.  Lastly,  let  Ri  be 
the  corresponding  element  of  the  diagonal  matrix  R /;  if  R /  were  not  diagonal,  a 
Cholesky  decomposition  [4:375]  would  need  to  be  computed.  This  decomposition 
would  re-define  the  filter  z  vector  and  the  H  matrix  but  would  then  give  a  diagonal 
R/  matrix  for  use  in  this  filter.  One  can  now  complete  the  update  by 


Then,  for  k  = 


f  = 

UrHr 

Vj  = 

Djj(  —  )fj  j  =  1,2, ...  ,n 

Oq  — 

R 

1,2,...  ,n 

0-k  = 

Ofc-l  +  fk^k 

Dkk(+)  = 

Dkk{-)dk-i/o-k 

Ok 

Vk 

Pk  = 

-fk/ak  -  1 

Ujk(+)  = 

Ujk{~)  +  bjPk  j  =  1,2, ...,(fc  —  1) 

bj  - 

bi  +  Ujk{-)vk  j  =  l,2,...,(fc  -  1) 

(2.63) 


(2.64) 
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The  filter  gain  is 


K  =  —  (2.65) 

a-n 

Once  the  gain  has  been  determined,  x(+)  and  P(+)  can  be  calculated  by  using 
Eqs  2.35  and  2.36. 

Upon  testing,  it  was  discovered  that  the  filter  covariance  would  not  converge; 
therefore,  the  filter  did  not  track  the  true  state.  The  maximum  eigenvalue  of  the 
covariance  and  the  true  error  were  growing  without  bound.  Further  investigation 
of  the  estimated  state  revealed  the  y  components  were  responsible  for  the  growth  in 
the  true  error.  As  the  state  was  propagated  through  several  orbits  the  downrange 
position  error  increased  to  more  than  1000  kilometers.  Therefore,  it  was  evident 
this  component  was  responsible  for  the  filter  instability. 

It  was  stated  earlier  in  this  section  that  the  U-D  factorization  filter  guaran¬ 
teed  stability.  This  statement  only  considered  numerical  stability,  it  said  nothing 
about  filter  stability.  In  order  to  guarantee  filter  stability,  the  system  must  be 
observable.  The  filter  performance  indicates  that  the  downrange  ( y )  component  of 
the  state  is  unobservable.  This  fact  seems  obvious  once  the  filter  has  revealed  it; 
the  downrange  displacement  of  the  satellite  cluster  from  the  reference  point  has  ab¬ 
solutely  no  effect  on  the  relative  range  measurements  or  on  the  relative  movement 
inside  the  cluster.  Mathematically,  this  is  verified  by  Eqs  2.45  through  eqxompzd. 
Remembering  the  note  on  the  initial  downrange  position,  it  can  be  seen  that  if  the 
entire  cluster  is  moved  downrange  from  the  reference  point  (by  a  constant  increase 
in  j/o  for  each  satellite),  neither  the  other  components  of  the  state  nor  the  range 
between  satellites  will  be  affected.  Since  the  measurements  contain  no  informa¬ 
tion  on  absolute  downrange  position,  it  is  impossible  to  estimate  the  absolute  y 
displacement  from  the  reference  point. 

2.3.3  Third  Filter.  This  instability  problem  can  be  resolved  with  the  U-D 
factorization  algorithm  but  the  filter  states  must  be  modified  in  order  to  eliminate 
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Figure  2.4.  Modified  State 

unobservability.  The  modified  state  replaces  the  absolute  downrange  component 
with  the  relative  downrange  displacement  Ay  as  illustrated  in  Figure  2.4.  Where 
Ay  is  defined  as 

Ay<  =  yi  -  yi  (2.66) 

For  t  =  2, 3, ....  s.  Therefore,  the  new  state  vectors  for  Xi  and  x,-  are: 


Because  the  states  have  changed,  the  $  matrix  must  also  be  modified.  By 
examining  Eqs  2.45  through  2.50,  one  sees  that  the  only  y,  and  therefore  the  only 
Ay,  dependance  is  in  Eq  2.46.  Substituting  Eq  2.46  into  Eq  2.66  yields 

Ay,(+)  =  [6(sin  xp  -  V»)](*i  -  *«)  +  Ay<  + 

2  ,  /  4  Zxp\ 

-  ( cos  xp  1 )  (i1-ii)+  -s  inV' - 1  (j/i  -  y.)  (2.68) 

IV  \V  V  ) 

where  the  state  components  on  the  right  side  of  the  equation  are  understood  to 

be  evaluated  at  the  previous  sample  time.  The  effect  of  this  modification  on  §i 

is  to  reduce  it  from  a  6-by-6  to  a  5-by-5  matrix  by  removing  the  row  and  column 

corresponding  to  the  y  component. 


4  —  3  cos  xp 

0 

3t7  sin  xp 
6r](cosrp  —  1) 


cos  ip 


—■qsinxp 


^(1  —  cost b)  0 


cos  ip  2  sin  xp  0 

— 2sinV’  —  3  +  4cosV>  0 

0  0  cos  ip 


(2.69) 


The  effect  on  3>,  is  that  the  second  row  of  the  matrix  is  now  the  negative  of  what 
it  was  in  Eq  2.51  The  second  element  (1)  being  the  exception  because  of  the  way 
that  Ay  was  defined  in  Eq  2.66. 


4  —  3  cos  xp  0 

— 6(sin  xp  —  xp)  1 

0  0 

3?7  sin  xp  0 

6rj(cosxp  —  1)  0 


-  z(cos  xp  -  1 )  -  4  sin  xp  +  0 


cos  xp 
—2  sin  xp 


2  sin  xp 
■3  +  4  cos  xp 


Unlike  the  structure  of  Eq  2.53,  when  the  individual  matrices  are  combined  into 
the  composite  state  transition  matrix,  the  result  is  not  a  block  diagonal  matrix. 
The  non  diagonal  elements  are 


6(s in^>  — V>)  0  jj(cosV>  —  1)  ^sinV’  —  ^  0 


= 


(2.71) 


Thus,  the  overall  system  is  of  the  form 


*(-)  = 


0  •••  0  xi(-) 

*ia  *9  •••  0  Xi(-) 


*t.  o  0  *,  I  I  x,(~) 


=  *x(+) 


(2.72) 


The  measurement  update  portions  of  the  filter  (h  and  H)  also  need  some 
revision.  The  non-linear  measurement  vector  h  requires  only  a  minor  alteration. 
Substituting  Eq  2.66  into  Eq  2.37  gives 


!{xx  -  x2)3  +  Ayi  +  (zi  -  z2)2 


'{x i  -  ®*)2  +  At/2  +  {zX  -  z,)2 


(2.73) 


The  linearization  of  this  vector  into  the  H  matrix  results  in  more  extreme  changes. 
Using  a  form  similar  to  Eq  2.39,  the  new  linearized  matrix  can  be  expressed  as 


Hx  Hx  0 
H2  0  H2 


H,-X  0  0 


(2.74) 


H.-X  j 
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where  for  t  =  1, 2, . . . ,  s  —  1 


0  0  0 


(2.75) 


*•  =  [  ~  (“t?11)  -t?  -(*1 H  0  0  °. 


(2.76) 


The  modification  is  complete  except  for  some  final  minor  changes. 

With  the  removal  of  the  second  diagonal  terms  of  the  Gj,  Po,  and  Qd  matri¬ 
ces,  U-D  covariance  factorization  algorithm  is  complete.  The  final  step  in  building 
this  filter  is  the  tuning  process  during  performance  evaluations. 
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III.  Performance  Analysis 


3.1  Analysis  Tools 

The  truth  model  and  the  U-D  factorization  covariance  filter  have  both  been 
developed.  All  the  necessary  information  for  an  initial  performance  analysis  is 
available;  specifically,  the  true  state  xt ,  t*~~  estimated  state  x ,  and  the  covariance 
P  have  all  been  computed.  The  remaining  task  is  to  convert  this  input  data  into 
useful  output. 

One  of  the  tools  for  this  data  conversion  is  the  true  error  et  as  depicted  in 
Figure  2.1.  The  quantities  of  interest,  y  and  yt,  can  be  determined  from  the  state 

by 

y*  =  ctxt  (3.1) 

and 

y  =  Cic  (3.2) 

Since  the  components  of  the  true  state  are  the  same  as  those  of  the  estimated  state, 

c  =  ct 

C  would  be  a  2-by-n  matrix  for  satellite  one  and  a  3-by-n  for  every  other  satellite 
(n  equals  the  total  number  of  states).  For  example, 


1  0  0  •••  0 

0  1  0  •••  0 

I 

.  Zt>  . 

.  . 

(3.3) 


and  the  true  error  is  given  by  Eq  2.1  as 


e«  =  yi  -  y*.  = 


(3.4) 


xi  —  xt, 

z\  -  *t, 


The  other  tool  needed  for  the  analysis  is  the  error  covariance.  If  this  system 
were  a  linear  system  driven  by  white  Gaussian  noise  with  zero  mean  processes, 
then  the  state  covariance  would  also  be  the  error  covariance  [4:page  335].  For  the 
system  in  this  study,  the  covariance  can  still  be  used  to  get  a  general  feel  for  the 
filter  performance  by  implementing  a  single  sample  Monte  Carlo  analysis.  For  a 
more  accurate  and  complete  analysis,  a  full  Monte  Carlo  simulation  would  need  to 
be  accomplished. 

In  this  scaled  down  study,  the  position  components  of  the  state  are  the  vari¬ 
ables  of  interest;  therefore,  the  corresponding  position  portion  of  the  error  covari¬ 
ance  Pe  is  needed  and  can  be  found  by 

=  CiPCf  (3.5) 

where  C<  is  the  same  matrix  used  in  Eqs  3.1  and  3.2,  and  where  i  designates  the 
satellite  of  interest. 

To  simplify  the  analysis  even  further,  the  magnitude  of  the  true  error  |eti  | 
is  plotted  against  the  standard  deviation  (o\)  of  the  covariance.  The  standard  de¬ 
viation  is  found  by  taking  the  square  root  of  the  largest  eigenvalue  of  Pe. .  These 
approximations  are  extremely  conservative  and,  therefore,  do  not  negate  the  valid¬ 
ity  of  the  performance  evaluation. 

3.2  Sample  Runs 

Sample  performance  runs  were  made  with  varying  numbers  of  satellites  in  the 
cluster  and  with  differing  values  of  Qj  and  P0.  All  filter  tuning  was  accomplished 
by  adjusting  Qj;  the  filter  measurement  covariance  matrix  R/  was  equal  to  Re  for 
all  time  and  was  based  on  assumed  range  measurement  errors  of  .01  meters.  The 
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Figure  3.1.  Filter  performance  (satellite  1  of  2,  no  dynamics  noise). 

cluster  radius  was  500  meters  at  a  reference  altitude  of  1000  kilometers  above  the 
earth.  This  resulted  in  a  cluster  orbital  period  of  approximately  6300  seconds. 

The  purpose  of  the  first  test  run  was  to  determine  if  the  filter  was  stable.  The 
system  was  initialized  to  two  satellites  with  the  estimator  knowing  true  position 
and  velocity.  The  dynamics  driving  noise  (Qj)  was  turned  off  and  the  filter  was 
fed  perfect  range  data.  Figure  3.1  shows  the  results.  As  hoped,  the  covariance 
converged  thereby  demonstrating  filter  stability.  The  divergence  of  the  true  error 
occurred  as  expected.  Since  the  filter  has  been  told  to  have  perfect  faith  in  its 
internal  dynamics  model  (Qd  =  0),  the  gain  will  converge  to  zero  and  decrease  the 
effect  of  the  range  measurement  inputs  on  the  updated  state.  Also,  note  the  true 
error  and  standard  deviation  are  periodic  with  the  orbit. 

Since  the  basic  stability  of  the  filter  was  confirmed,  the  next  test  was  to  check 
for  filter  performance  improvements  with  the  addition  of  dynamics  driving  noise. 
The  order  of  magnitude  of  the  Qj  matrix  diagonals  was  determined  by  calculating 
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Figure  3.2.  Filter  performance  (satellite  1  of  2,  with  dynamics  noise). 

the  true  error  after  the  initial  time  propagation.  The  values  used  were  10-11  for 
position  and  10-17  for  velocity.  The  resulting  plot  (Figure  3.2)  reflects  a  stable,  well 
behaved  filter.  However,  the  values  used  were  apparently  too  large.  The  desired 
end  result  of  tuning  is  to  have  agreement  between  the  filter’s  internally  generated 
error  estimates  said  the  true  error.  In  this  case,  the  filter  is  overestimating  the 
errors  it  is  committing  and  putting  too  much  weight  on  the  measurements  [4:page 
338].  Another  problem  is  the  magnitude  of  the  true  error;  two  kilometers  is  well 
outside  the  allowable  error  tolerances  outlined  in  Appendix  A.  Further  tuning  is 
required. 

Another  method  for  decreasing  the  true  error  is  by  increasing  the  cluster 
size.  Each  additional  satellite  provides  more  information  to  the  state  estimator 
consequently  improving  the  position  estimate  of  satellite  one.  Figure  3.3  shows 
the  dramatic  results  of  increasing  the  array  to  five  elements.  The  initial  transient 
is  much  smaller  as  are  the  «r  and  e(  values.  However,  the  true  performance  is 
revealed  when  the  scale  size  is  expanded.  Figure  3.4  depicts  a  filter  significantly 
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underestimating  its  own  erro~'  to  the  point  of  ignoring  the  measurement  data. 
Even  though  more  diagonal  ter  -  were  added  to  the  dynamics  noise  covariance 
with  each  additional  satellite,  it  still  was  not  enough  to  keep  the  filter  gain  from 
going  to  tero.  The  probable  cause  of  this  effect  is  significant  unmodeled  cross¬ 
correlation  effects  in  the  Q4  matrix.  The  time  scale  on  the  latest  figure  begins  at 
3500  seconds  because  the  transient  respr  ise  of  the  filter  was  off  the  scale  of  the 
plot.  Since  the  filter  appears  to  be  in  steady-state  within  three  orbital  periods, 
future  plots  will  terminate  at  18900  seconds. 

The  objectives  for  the  next  filter  test  were:  to  reduce  the  tran^i-nt  response; 
to  increase  the  cluster  to  its  maximum  sise  of  ten  satellites;  and  to  do  a  rough 
retuning  in  order  to  evaluate  performance.  The  transient  response  was  improved 
by  replacing  the  diagonal  elements  of  the  Po  matrix  with  those  of  the  steady-state 
covariance  matrix.  The  new  diagonals  were  on  the  order  of  10~*  for  position  and 
10-1’  for  velocity.  Also,  the  dynamics  noise  covariance  was  increased  to  2  X  10“* 
and  2  x  10~M  for  position  and  velocity  variances  respectively.  Figures  3.5  and  3.6 
show  the  results  of  these  changes.  The  ten  satellite  filter  perforir  ed  exceptionally 
well.  The  transient  response  of  the  standard  deviation  was  reduced  from  a  peak 
of  approximately  3  kilometers  to  a  maximum  of  around  2.5  meters  (Figure  3.6). 
The  maximum  true  error  was  one  tenth  of  a  meter  for  satellite  one  and  close  to 
zero  for  the  others.  It  seemed  reasonable  to  expect  satellite  one’s  position  to  be 
known  best  because  all  the  range  measurements  to  every  cluster  element  included 
information  on  this  single  satellite.  The  only  apparent  explanation  lies  in  the 
differences  between  the  state  model  of  satellite  one  and  the  model  of  the  other 
satellite  states  as  seen  in  Eq  2.67. 

On  this  run,  all  the  filter  statistics  remain  the  same  as  in  the  previous  test. 
The  only  change  is  that  now  corrupted  data  will  be  fed  into  the  estimator.  In 
the  previous  case,  with  perfect  measurements,  the  true  error  was  less  than  one 
tenth  of  a  meter.  The  true  error  response  to  range  data  corrupted  at  .01m  is 
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Figure  3.5,  Filter  performance  (satellite  1  of  10,  with  dynamics  noise  and  steady 
state  covariance). 


Figure  3.6.  Filter  performance  (satellite  10  of  10,  with  dynamics  noise  and  steady 
state  covariance). 


Figure  3.7.  Filter  performance  (satellite  1  of  10,  with  dynamics  noise,  steady- state 
covariance,  and  noisy  data). 

seen  in  Figures  3.7,  3.8,  and  3.9.  The  results  were  in  agreement  with  what  was 
anticipated.  The  true  errors  grew  larger.  In  the  case  of  satellite  number  one,  the 
error  grew  to  the  point  that  it  exceeded  the  internally  estimated  error.  But  this 
could  be  corrected  by  again  retuning  the  filter.  Even  without  retuning,  the  total 
error  was  still  less  than  one  meter  which  is  well  within  the  minimum  required 
accuracy. 

Finally,  as  a  follow  on  to  the  previous  test  run,  the  filter  was  retuned  to 
improve  performance.  The  dynamics  noise  covariance  was  increased  to  8  x  10~* 
and  8  x  10H  for  position  and  velocity  variances  respectively.  Figures  3.10,  3.11, 
and  3.12  show  that  the  performance  was  improved  by  minimal  tuning,  as  claimed. 


Figure  3.10.  Filter  performance  (satellite  1  of  10,  retuned,  with  noisy  data) 
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Figure  3.11.  Filter  performance  (satellite  5  of  10,  retuned,  with  noisy  data) 


Figure  3.12.  Filter  performance  (satellite  10  of  10,  retuned,  with  noisy  data) 


IV.  Conclusion 


The  U-D  factorization  filter  performed  well  for  all  cases  evaluated.  Even 
though  tuning  was  limited  to  the  dynamics  noise  matrix,  the  estimator  was  able  to 
significantly  exceed  the  minimum  required  accuracy.  Numerical  stability  and  filter 
stability  /ere  maintained  during  all  runs  and  the  initial  transients  were  successfully 
reduced  by  many  orders  of  magnitude.  All  of  this  was  accomplished  with  a  cluster 
size  significantly  larger  than  previously  tested.  The  analysis  performed  in  this 
thesis  was  only  an  initial  evaluation  there  remains  much  still  to  do. 

A  full  Monte  Carlo  analysis  could  be  conducted  to  determine  the  actual 
filter  error  statistics  including  biases.  Minimize  these  errors  through  tuning.  In 
the  tuning  process,  investigate  the  performance  enhancement  in  adjusting  non¬ 
diagonal  terms  in  the  Q *  matrix.  Do  not  limit  the  R/  matrix  to  the  values  of  the 
truth  model. 

Test  the  robustness  of  the  filter.  Analyze  how  accurate  the  initial  conditions 
on  the  state  must  be  for  the  estimator  to  lock  on.  Test  how  noisy  the  range  data 
can  be  and  still  obtain  valid  estimates.  How  far  can  the  cluster  be  dispersed  before 
the  assumptions  in  the  Clohessy- Wiltshire  equations  become  invalid  and  nullify  the 
mathematical  basis  of  the  filter?  Finally,  test  the  filter  against  a  more  accurate, 
real  word  truth  model. 
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Appendix  A.  Satellite  Precision  Requirement 


A.l  Introduction 

This  appendix  outlines  how  the  position  accuracy  requirements  were  deter¬ 
mined.  These  requirements  were  used  to  evaluate  filter  performance. 


A.S  Accuracy  Requirement 

As  stated  in  the  main  introduction,  the  positions  of  the  satellites  must  be 
known  to  at  least  one  quarter  of  the  radar  wavelength.  The  maximum  wavelength 
is  a  function  of  the  size  of  the  antenna. 

According  to  the  problem  description,  the  cluster  consists  of  ten  satellites  at 
an  altitude  of  1000  km.  Assume  that  each  satellite  has  a  ground  coverage  area 
approximating  that  of  Wyoming.  Estimate  the  radius  of  this  area  as  480  km. 
Therefore,  the  angle  subtended  at  the  satellite  by  the  ground  is 


0- 


defining  this  angle  as  the  3-dB  beamwidth  and  solving  for  it  gives 


03,b  =  2  tan-1  (^) 


(A.l) 


(A.2) 


Assuming  that  the  antenna  efficiency  £  is  55  percent,  the  antenna  gain  is 
found  by  [7:page  81] 

G  “  <A-3) 


If  the  conservative  estimate  is  made  that  the  overall  gain  is  just  the  multiple  of  the 
individual  gains  then 

Gtot  =  10  x  G  (A. 4) 


The  radar  wavelength  can  then  be  found  by 

Air  A 

Lttot  =  q— rr~ 


(A.S) 


where  A  is  the  cluster  area.  Finally,  solving  for  A 

A  =  ,r(-i-r  (A.6, 

\LttOT  } 

where  r  is  the  cluster  radius.  Remember,  the  needed  accuracy  is  at  least  twenty-five 
percent  of  this  wavelength.  Table  A.l  gives  some  representative  values  of  minimum 
required  filter  performance. 


cluster  radius 

A/4 

.5  km 

5  km 

50  km 

2.75  m 
275  m 
2750  m 

Table  A.l.  Cluster  radius  vs  required  accuracy 
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